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An accurate and efficient solution of the radiative transport equation is proposed for modeling the 
propagation of photons in the three-dimensional anisotropically scattering half-space medium. The exact 
refractive index mismatched boundary condition is considered and arbitrary rotationally invariant 
scattering functions can be applied. The obtained equations are verified with Monte Carlo simulations in the 
steady-state, temporal frequency, and time domains resulting in an excellent agreement. 

I n many research fields the radiative transport equation (RTE) is the fundamental equation for describing 
I particle transport like neutron propagation in reactor physics 1 or photon propagation in atmospheric physics, 
I in astronomy, or in biophotonics 2 " 4 . For quantitative investigations of the photon transport in biological media 
the RTE is used for all scales, i.e. for macroscopic, mesoscopic, and microscopic applications 5 , where interference 
effects don't have to be taken into account. For example, it is applied for non-invasive investigations of the light 
propagation in the human brain, in small animals, or for quantitative microscopy 6-10 . 

The RTE is usually solved with numerical methods. In the case of isotropic scattering half-space media the RTE 
can be solved analytically as it is reported in the publications 11 " 14 . In the field of biophotonics mainly the Monte 
Carlo (MC) method is applied 15 . Due to the long simulation time needed with numerical methods until accurate 
results are received the RTE is often approximated by the diffusion equation. Within the diffusion approximation 
(DA) analytical solutions can be obtained for several geometries 16 " 19 . However, it is well-known that the DA fails 
in the high spatial and temporal frequency domain (TFD) resulting in a breakdown of the diffusion theory for 
small distances to sources and boundaries as well as for short time values. The fact that adequate analytical 
solutions are still rare but of great interest is underlined by many recently published articles which present 
diffusion theory based models 20 or approximated solutions of the RTE 21 . 

In this article we report on radiative transport solutions for describing the propagation of photons in the three- 
dimensional anisotropically scattering half-space medium. The effect of Fresnel reflection at the interface which is 
present in nearly all measurements is also considered. The semi- infinite geometry is abundantly used in all three 
measurement domains 16,17-22,23 . The accuracy of the derived equations is demonstrated in the figures below by 
comparisons with Monte Carlo simulations showing an excellent agreement. 

Results 

In the following we give an overview regarding the applied analytical approaches and the appropriate numerical 
calculation steps for obtaining the proposed radiative transport solution. We start with the derivation in the TFD 
where the quantities of interest are the demodulation and the temporal phase shift between the injected specific 
intensity and the measured signal. 

The specific intensity caused by a sinusoidally modulated pencil beam 

Imcip, 8 t) = S[p)5[a-i)ap[icot) (1) 

which enters the boundary of the semi-infinite scattering medium z s 0 becomes G(r, s, t) = G(r,s) exp(tor), 
where p e R 2 and m e R is the angular modulation frequency. The complex amplitude obeys the RTE 2,3 



s-VG(r, s)+n t G(t, s)=pi s l G(r, s')/(s-s')ds' 



(2) 



where \i t = [i a + [i s + ico/c is the frequency dependent total attenuation coefficient, [i a the absorption coefficient, 
[i s the scattering coefficient, c = c 0 ln is the velocity of light in the scattering half-space with refractive index n. The 
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detection of the specific intensity takes place at the spatial position 
reR 3 , whereas the unit vector se§ 2 specifies the direction of the 
photon propagation and /(s - s') is the scattering phase function 
which is normalized to unity. The exact boundary condition (BC) 
corresponding to a pencil beam which enters the half-space z s 0 is 
defined in the angular domain s-z > 0 and is given by 

G(p, 0, fi, </>)=R(, J .)G{p, 0,-fi, <j>) + d(p)8(s-z), (3) 

where R(fi) is the Fresnel reflection coefficient 12 and /( = cosO. Here 
Green's function in the frequency domain is separated into its bal- 
listic and diffuse contribution 3 according to G(r,s) = G(r,s) + /(r,s). 
By substituting this ansatz in (2) we obtain the first RTE 

s-VG fc (r,s)=-/i t Gi,(r,s), (4) 

whose solution under consideration of the BC (3) yields the ballistic 
component 

G fc (r, s)=8(p)e-''- z e(z)S(s-i), (5) 

where 0(z) denotes the Heaviside step function. Taking into account 
this component we obtain the second RTE for the diffuse part 



s-VI(r, s) +^/(r, g) = n s f I(r, S')/(S-s')dS' 

Js> 2 

+ ft <5(/>) e -"< z ©( Z )/(s-z), 



which has to be solved subject to the BC 

I(p, 0, fi, <j})=R(n)I(p, 0,-fi, < 



s-z>0. 



(6) 



(?) 



The particular solution of (6) can be evaluated rapidly by performing 
a similarity ansatz. It can be shown that 

lW(r, s)=e-"- z ®(z)J2lUp)Y, m (0, x) (8) 



with Hankel transformed coefficients 



4n0)= — 
In 



1 lin {q)Jm{qp)qdq 



(9) 



is a complete and exact solution of the RTE (6) for appropriate 
spectral components Ii m (q) = {~l) m h,-m{q)- Moreover Yi m (0, j) 
with x = <t> ~ §p are the spherical harmonics (SH) and J m (x) are 
the Bessel functions of the first kind. Inserting (8) together with (9) in 
(6) and making use of recurrence relations for the SH and the Bessel 
functions leads to the following set of linear equations labeled by Z = 
m,m + 1, ... for m = 0, 1, ... according to 



-/'/ 



V2T- 



•>(<?) + 



M) 



(10) 



■1 ' ~" " \ 2/ • 3 
+ iq[bi- i,mh-i,m+i (q) - Cl+l, m Il+ i, m+ 1 [q)] 
+iq[h + l,m-*/+ l,m — 1 (q) — c l-l.mh-l,m-l 

(«)] 

+ V2l+la,I, m (q) = (2/+ 1)/,<W S / v 7 ^, 



with the quantities cr ; = fi t — fifi s . The remaining quantities are given by 
ai m = Vl 2 -m 2 , b lm = c t _ m and c, m = ^/(l + m)(l + m+ l)/(2/+ 1)/ 
2. The expansion coefficients of the scattering phase function are 
defined as 



fl=l /( 
Js 2 



s-s')P/(s-s')ds' 



(11) 



where Pi(x) are the Legendre polynomials. The system of linear equa- 
tions (10) written in matrix notation becomes 

(-H t A+iqB+S)\I(q)) = \e). (12) 

Now the solution to the boundary-value problem (7) can be achieved 
by considering additionally the general solution to the homogeneous 
RTE. In our recent study 26 we solved the homogeneous problem in the 



steady-state domain by applying the method of rotated reference 
frames 24,25 . Although the obtained solution is formally exact its success- 
ful application is impossible for several situations of high practical 
importance due to numerical shortcomings. Therefore, we reconsid- 
ered the homogeneous problem in the TFD and obtained the following 
complete new result 



/W(r 



$> im (p, z)Y lm (0, x), 



(13) 



where 



with 



^im(P^ z } = J n J o ^im^ z)Jm{qp)qdq (14) 



"Alm^ z)= C '("3)exp(-C,z)^(q) 



(15) 



a,->0 



and ctj = Re(c,). The eigenvalues c, and the eigenvector components 
^hl (?) f°ll° w as solution of the generalized eigenvalue problem 



(S+iqB)\u(q)) = SA\u(q)), 



(16) 



which can be translated into a standard eigenvalue problem because 
det(S + iqB) 0. The characteristic equation of (16) is given by 



n ri-^)=o, 



(17) 



where A ( are the eigenvalues of the matrix (S + iqBY 1 A. Note that in 
the steady-state domain which entails that co = 0 the eigenvalues C ; = 
II A t are strictly real- valued. The corresponding nontrivial solution 
I u i(q)) ^ 0 to the eigenvalue ^ is that which is given by the eigenvector 
to X t , The general solution to the RTE (6) for the semi-infinite geometry 
is now obtained via superposition of the homogenous and particular 
part. 

Next, we apply the generalized Marshak BC 24,25,27 for satisfying 
equation (7) and to determine the unknown constants C t {q) within 
the homogenous solution. At this stage the solution to the RTE (6) for 
the semi-infinite geometry is completed. In view of diagnostic and 
therapeutical applications the back-scattered light R(p) = — f f-d{/>, 0, 
s)ds and the internal fluence <l>(r) = / J(r, s)ds are quantities of 
particular interest. By making use of the derived specific intensity 
as well as of the orthogonality relation jYi m (s)Y^ m ,(s)ds = <5/f(5 mm ', 
we find the following expressions for the reflectance 



R(p) = 

and for the fluence in z s 0 



3 (p,0)+/ 10 (p)] (18) 



<D(r) = 5(p)e-^ + V4^# 00 (p, z) + e -"' z /oo(p)] • (19) 

The above derived solution can be evaluated for an arbitrary angular 
frequency. Therefore, we can reconstruct the time-resolved impulse 
response or Green's function for the diffuse specific intensity accord- 
ing to 



I[t, s, r)= — ( I(r, s, oj)e ioJ 'dco. 
2nj 



(20) 



The frequency dependent coefficients from system (10) satisfy the 
condition Ii m (q, w) =IT_ m (q, —co). In addition, the real and ima- 
ginary parts of J(r, s, co) are connected by the Hilbert transform so 
that it is sufficient to invert e.g. only the real part. Green's function of 
the time-dependent RTE for the semi-infinite medium is com- 
plete by adding the ballistic contribution Gf,(r, s, t) = ce~''' ct ©(r) 
<5(r — dz)(5(s — z). 

The derived solutions to the RTE are verified by comparisons with 
MC simulations 28 in three different measurement domains. The 
lateral intensity profile of the collimated incident beam is modeled 
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p (mm) 

Figure 1 | Steady- state reflectance versus radial distance p due to an incident unmodulated light source having two different beam radii if. The 

anisotropy factor of the HG phase function is g = 0.8. 



by a Gaussian function according to S(p) = 2/{mf) exp( — 2p 2 /if), 
where denotes the radius of the beam. Due to the convolution 
theorem we have to multiply the right-hand side of system (12) with 
the Hankel transformed distribution S(q) = exp( — q 2 if/8). The bal- 
listic specific intensity becomes now Gf,(r, s) = S(p)e~' ,|Z 0(z) 
<5(s — z). In most cases shown below the absorption coefficient and 
the reduced scattering coefficient of the semi-infinite medium are 
assumed to be jx a = 0.01 mm -1 and |/j=(l-j)fi s =l.(l mm -1 , 
respectively, which are typical values for biological tissue in the 
near-IR spectral range. If not stated otherwise we consider the 
Henyey-Greenstein (HG) scattering phase function 4 for different 
values of the anisotropy factory, where the corresponding expansion 
coefficients are given by/; = g 1 . The relative refractive index of the 
non-scattering medium is set to 1.0 whereas the refractive index of 
the scattering medium is set to n = 1.4 if not stated otherwise. In the 
following we refer to the developed transport solution as the modi- 
fied spherical harmonics method (MSHM). 

We start with the spatially-resolved reflectance in the steady-state 
domain which is caused by an incident unmodulated light source 
having two different beam radii }]. Figure 1 contains the reflectance 
obtained from equation (18) as well as the result of the MC method. 



An excellent agreement is obtained for both beam radii. Thus, we also 
computed the relative differences between both methods in the case 
of the smaller beam radius and included them in the inset. Within the 
stochastic nature of the MC simulations an exact agreement is found. 

Figure 2 displays the steady-state reflectance as function of the 
radial distance for a semi-infinite medium having two different 
absorption coefficients. For the less absorbing tissue we computed 
the relative differences between the MC simulation and the MSHM 
which are depicted in the inset. In accordance with the first compar- 
ison the MSHM leads practically to the same reflectance as the MC 
simulations for both absorption coefficients. 

For the next comparison the absorption coefficient of the scatter- 
ing half-space is significantly increased relative to the figures above. 
Figure 3 shows the steady-state reflectance obtained from the MC 
simulation (symbols), the MSHM (solid lines) and the DA (dashed 
lines). In that cases the DA cannot be applied for modeling ade- 
quately the back- reflected photons which is a well-known shortcom- 
ing of the diffusion theory. 

The influence of the scattering phase function on the spatially 
resolved reflectance is considered in figure 4. To this end we use 
the HG phase function with g = 0.9 for the case of a highly forward 
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Figure 3 | Spatially-resolved reflectance as function of p due to an 
incident unmodulated beam with radius 1/ = 0.3 mm evaluated for two 
different absorption coefficients. The anisotropy factor of the HG phase 
function is g = 0.8. 
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Figure 4 | Spatially-resolved reflectance as function of p due to an 
incident unmodulated beam with radius tj = 0.3 mm evaluated for two 
different scattering phase functions. The anisotropy factor of the HG 
phase function is g = 0.9. 



scattering medium and the Rayleigh function /(s - s') = 3[l + (s - 
s') 2 ]/(167i) for the case of almost isotropic scattering. Within the 
MSHM this function is given using the expansion coefficients// = 
(5/o + 0.1(5/2. The solid lines correspond to the MSHM whereas the 
symbols display the results of the MC simulations. In addition, the 
reflectance in DA 29 caused by the assumed Gaussian beam is repre- 
sented by the dashed line. Again, the agreement between the analyt- 
ical and the numerical solution of the RTE is excellent, whereas DA 
performs bad especially at short distances. 

We now consider the Mie phase function i.e. the solution of 
Maxwell's equations for the scattering of a plane wave by a sphere 30 . 
To this end we computed the expansion coefficients according to Eq. 
(11). In figure 5 the spatially-resolved reflectance for the Mie phase 
function and the Rayleigh function are shown and compared to the 
results obtained by the MC method (symbols). For calculation of the 
Mie phase function a sphere having a diameter of 1 fim and a refract- 
ive index of 1.59 is assumed to be surrounded by water with index 
1.33. The wave length of the unpolarized incident plane wave is 
assumed to be k = 633 nm. The refractive index of the scattering 
medium is set to n = 1.33. 

Figure 6 displays the internal steady-state fluence according to 
equation (19) due to an unmodulated narrow incident Gaussian 



r 




p (mm) 



Figure 6 | Steady-state fluence versus radial distance p due to an incident 
unmodulated light source. The anisotropy factor of the HG phase 
function is g = 0.8. 




p (mm) 

Figure 5 | Spatially-resolved reflectance as function of p due to an incident unmodulated beam with radius // = 0.8 mm evaluated for the Mie phase 
function and the Rayleigh function. 
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Figure 7 | Reflectance amplitude and temporal phase shift versus radial distance p due to an incident modulated beam with radius tj = 0.3 mm. The 

anisotropy factor of the HG phase function is g = 0.6. 



beam with radius r\ = 0.05 mm. The fluence is evaluated very close to 
the reflecting boundary at z = 0.05 mm. Similar as in the compar- 
isons above the MSHM (solid line) coincidences excellently with the 
internal fluence obtained from the MC simulation (noisy curve and 
open circles in the inset). 

In figure 7 we show a comparison between the MSHM and the MC 
method in the TFD. To this end we evaluated the reflectance R(p) = 
A(p) exp[ — i<p(p)] according to equation (18) for two values of the 
angular modulation frequency w = 2nf which covers the typical 
frequencies used in diffuse optical imaging. The amplitude A(p) = 



\R(p)\ and the phase shift <p(p) = — arg[i?(p)] are verified with MC 
simulations. Here the corresponding MC results are obtained via a 
discrete Fourier transformation of the time-resolved reflectance 
caused by an infinitely short light pulse. The inset shows the phase 
shift for the smaller source-detector separations which exhibits an 
noticeable curved course due to the Gaussian intensity profile of the 
incident beam. It can be seen that the MSHM enables an accurate 
modeling of the sinusoidally modulated reflected light. 

For the next comparison we evaluated equation (18) for a variety 
of angular modulation frequencies to reconstruct the time-resolved 



OH 




0.2 0.3 
t(ns) 



0.5 



Figure 8 | Time-resolved reflectance due to an infinitely short light pulse evaluated for two different source-detector separations p. The anisotropy 
factor of the HG phase function is g = 0.6. 
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impulse response for the reflected light R(p, t) according to equation 
(20). Figure 8 shows the time-resolved reflectance for two different 
source-detector separations. Here the solid lines correspond to the 
MSHM, the noisy curves (and the symbols in the inset) are the result 
of the MC simulations and the dashed lines denote the solution of the 
diffusion equation under the extrapolated BC 4 . The radius of the 
lateral intensity profile is f] = 0.5 mm. While the MSHM and MC 
simulation are in excellent agreement for both source-detector 
separations the DA shows significant differences compared to the 
transport theory solutions even for relatively large time values. 

For a fixed source-detector separation of p = 4.05 mm we inves- 
tigated the influence of the refractive index mismatch on the time- 
resolved reflectance. Figure 9 displays the result of the MSHM (solid 
lines), the MC simulation (noisy curves, symbols in the inset) and the 
DA (dashed lines). The radius of the lateral intensity profile is f\ = 
0.3 mm. 

Similar to the previous comparison the DA cannot be applied for 
accurate modeling the photon propagation for both matched and 
mismatched BC for early times. 

Discussion 

In conclusion, we have proposed a radiative transport solution for 
modeling the photon propagation in the three-dimensional semi- 
infinite medium which is an essential advance in view of the 
currently available solutions and approaches for solving the three- 
dimensional RTE. To this end the RTE was solved subject to the exact 
specular reflecting BC. The particular solution to the inhomogeneous 
RTE was obtained without recourse to the spatially- resolved impulse 
response by seeking a similarity solution which can be evaluated 
rapidly and accurately. In contrast to the classic spherical harmonics 
method (P N approximation) the proposed solution exhibits decisive 
differences and novelties. For example, the system of linear equations 
(12) as well as the eigenvalue problem (16) depend only on the scalar 
wave number q instead of the two-dimensional vector q = (q lt q 2 ) 
resulting in considerable simplifications and reduced computational 
cost. In addition, the number of equations is reduced by a factor of 2. 
The derived equations were successfully verified by comparisons 
with MC simulations in the three commonly measurement domains 
showing, within the stochastic nature of the simulations, an exact 
agreement. 

Finally, we mention that the proposed method can be extended 
in several ways. For example, arbitrary spatial and temporal shapes 
of the incident source can be considered and the solid angle for 
calculating the reflectance can be changed. A further important 



extension is that solutions for more complicated geometries such 
as the layered tissue structure can be obtained without alteration of 
the overall solution approach which will be presented in a forthcom- 
ing article. 

Methods 

Implementation and evaluation. The proposed solution was implemented with a 
non-optimized MATLAB script. The zero order Hankel transform which appears in 
Eqs. (18) and (19) was evaluated using the Gaussian Legendre quadrature method 
with 240 supporting points. The time-resolved reflectance was generated from the 
frequency domain Green's function via the discrete Fourier transform. The 
computation time needed for obtaining the results in the steady-state and the TFD 
was about one second whereas the time-resolved reflectance required about 30 
seconds. The corresponding MC simulation time took about 20 hours for the steady- 
state domain and about four days in the TFD and time domains. 

We additionally note that due to convenience the specific intensity of the incident 
beam which has entered the scattering medium was normalized to unity. This means 
that the specular reflection at the entry into the semi-infinite medium has not been 
considered. However, this effect can be easily taken into account by multiplying the 
final results above with the transmission coefficient T — 1 — (n — l) 2 /(n + l) 2 . 

The ballistic contribution. In order to derive the ballistic contribution we can 
consider the inhomogeneous version of Eq. (4) 

s-VGi,(r, s) + /i,G fc (r, s) =<5(r)<5(s-z). (21) 

One possibility to solve this equation is given by performing the Fourier transform 
Gf,(k, s) = /Gf,(r, s)e~' kr d 3 r. Therefore, the spectral component for the infinite 
medium becomes 



whose inverse is a standard transform pair which corresponds with Eq. (5). It can be 
seen that the infinite-space Green's function is simultaneously the wanted half-space 
Green's function which satisfies the specular reflection BC (3). 
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